#####
# Script to generate figure 1. Labeling was done after image generation using Google Slides.
#####

p <- 1
r_s <- 0.05
r_l <- 0.05
F <- 0.35
bar_S <- 5
b <- 1/bar_S
g_0 <- 1.25 # this is sigma in the paper. For scenario where only open-access causes Kessler
g_1 <- 20 # For scenario where both open-access and planner cause Kessler

# Survival probability function
survival <- function(S,b){
	pmax(1 - b*S,0)
}

# Derivative of survival probability function with respect to S
survival_prime <- function(S,b){
    (-b)*(S < 1/b) + 0*(S >= 1/b)
}

# Long-run object stock
X_l <- function(S,b,g) {
	value = S + g*(1-survival(S,b))*S
	return(as.numeric(value))
}

# Derivative of long-run object stock with respect to S
X_l_prime <- function(S,b,g) {
    value = 1 + g*(1-survival(S,b)) - g*b*S
    return(as.numeric(value))
}

# Expected marginal revenues
expected_MR <- function(S,p=p,r_s=r_s,r_l=r_l,b=b,g=g){
	p*(survival(S,b)/(1+r_s))*( 1 + survival(X_l(S,b,g),b)/(1+r_l) )
}

# external cost
EC <- function(S, p, r_s, r_l, b, g) {
    -(p/(1+r_s)) * ( survival_prime(S,b) * (1 + survival(X_l(S,b,g),b) / (1+r_l)) + survival(S,b) * ( survival_prime(X_l(S,b,g),b) * X_l_prime(S,b,g)) / (1+r_l) )
}

# Kessler threshold
S_K <- function(bar_S, g) {
    obj  = function(S,g,bar_S) {X_l(S,1/bar_S,g) - bar_S}
    uniroot(
        f = obj, g=g, bar_S=bar_S, interval=c(0,bar_S)
    )$root
}

##### Function test blocks
## Scenario 0: Only open-access causes Kessler
# Calculate open-access equilibrium solution. Uncomment to play around with uniroot to explicitly compute the solution.
# uniroot(
#     f = function(S, p, r_s, r_l, b, g_0, F) {
#         expected_MR(S,p,r_s,r_l,b,g_0) - F
#     }, p=p, r_s=r_s, r_l=r_l, b=b, g=g_0, F=F, interval=c(0,bar_S)
# )$root

# Calculate planner solution. Uncomment to play around with uniroot to explicitly compute the solution.
# uniroot(
#     f = function(S, p, r_s, r_l, b, g_0, F) {
#         expected_MR(S,p,r_s,r_l,b,g_0) + EC(S,p,r_s,r_l,b,g_0) - F
#     }, p=p, r_s=r_s, r_l=r_l, b=b, g=g_0, F=F, interval=c(0,bar_S)
# )$root

# Calculate Kessler threshold
# S_K(bar_S,g_0)

# ## Scenario 1: Both open access and planner cause Kessler
# # Calculate open-access equilibrium solution. Uncomment to play around with uniroot to explicitly compute the solution.
# uniroot(
#     f = function(S, p, r_s, r_l, b, g_1, F) {
#         expected_MR(S,p,r_s,r_l,b,g_1) - F
#     }, p=p, r_s=r_s, r_l=r_l, b=b, g=g_1, F=F, interval=c(0,bar_S)
# )$root

# # Calculate planner solution. Uncomment to play around with uniroot to explicitly compute the solution.
# uniroot(
#     f = function(S, p, r_s, r_l, b, g_1, F) {
#         expected_MR(S,p,r_s,r_l,b,g_1) + EC(S,p,r_s,r_l,b,g_1) - F
#     }, p=p, r_s=r_s, r_l=r_l, b=b, g=g_1, F=F, interval=c(0,bar_S)
# )$root

# # Calculate Kessler threshold
# S_K(bar_S,g_1)
#####

# Calculate Kessler thresholds
S_K_0 <- S_K(bar_S,g_0)
S_K_1 <- S_K(bar_S,g_1)

# Run script
lapply(c("tidyverse", "ggplot2", "patchwork", "ggthemes"), library, character.only=TRUE)

# Create a grid of S values and associated launch rates
S_grid <- data.frame(S=seq(0,bar_S,length.out=1000)) %>%
    mutate(
        X_l_0 = X_l(S,b,g_0),
        X_l_1 = X_l(S,b,g_1)
    )

# Generate values for scenario 0
output_g0 <- S_grid %>%
	mutate(EMR = expected_MR(S,p=p,r_s=r_s,r_l=r_l,b=b,g=g_0),
		PMC = F,
        EC = EC(S,p=p,r_s=r_s,r_l=r_l,b=b,g=g_0),
		SMC = PMC + EC(S,p=p,r_s=r_s,r_l=r_l,b=b,g=g_0),
		surviving_sats = S*survival(S,b=b),
		ETR = S*EMR,
		PTC = PMC*S,
		STC = SMC*S,
        S_K = S_K_0
		) 

# Generate values for scenario 1 
output_g1 <- S_grid %>%
	mutate(EMR = expected_MR(S,p=p,r_s=r_s,r_l=r_l,b=b,g=g_1),
		PMC = F,
        EC = EC(S,p=p,r_s=r_s,r_l=r_l,b=b,g=g_1),
		SMC = PMC + EC(S,p=p,r_s=r_s,r_l=r_l,b=b,g=g_1),
		surviving_sats = S*survival(S,b=b),
		ETR = S*EMR,
		PTC = PMC*S,
		STC = SMC*S,
        S_K = S_K_1
		) 

# Join the two outputs 
output <- left_join(output_g0, output_g1, by="S", suffix=c("_g0","_g1"))

# Plot with S on the x-axis, and EMR, PMC, SMC on the y-axis. There is a dotted vertical line showing the Kessler threshold. The EMR is in solid black, the PMC is in solid black, the SMC is in solid red, and the Kessler threshold is in dashed black. The plot is on the left, and the legend is on the right. Plot is made with ggplot. There are two panels. On the left, with output_g0. On the right, right output_g1.
## Scenario 0 plot
simple_model_plot_g0 <- ggplot(data=output_g0, aes(x=S)) +
    geom_line(aes(y=EMR), color="black") +
    geom_line(aes(y=PMC), color="black") +
    geom_line(aes(y=SMC), color="blue") +
    geom_vline(xintercept = S_K(bar_S,g=g_0), linetype="dotted") +
    labs(x = "S", y = "Revenues and costs", title = "") +
    ylim(c(0,max(output_g0$EMR))) +
    scale_x_continuous(name = "S", breaks = c(bar_S), label = expression(bar(X))) +
    theme_minimal() + 
    theme(
        axis.text.y = element_blank(),
        axis.text.x = element_text(color='black'),
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
        ) +
    geom_segment(x=0,y=0, xend=5,yend=0, size=0.15) + geom_segment(x=0,y=0, xend=0,yend=5, size=0.15) # axis bars
## Scenario 1 plot
simple_model_plot_g1 <- ggplot(data=output_g1, aes(x=S)) +
    geom_line(aes(y=EMR), color="black") +
    geom_line(aes(y=PMC), color="black") +
    geom_line(aes(y=SMC), color="blue") +
    geom_vline(xintercept = S_K(bar_S,g=g_1), linetype="dotted") +
    labs(x = "S", y = "Revenues and costs", title = "") +
    ylim(c(0,max(output_g1$EMR))) +
    scale_x_continuous(name = "S", breaks = c(bar_S), label = expression(bar(X))) +
    theme_minimal() + 
    theme(
        axis.text.y = element_blank(),
        axis.text.x = element_text(color='black'),
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
        ) +
    geom_segment(x=0,y=0, xend=5,yend=0, size=0.15) + geom_segment(x=0,y=0, xend=0,yend=5, size=0.15) # axis bars

simple_model_plot <- simple_model_plot_g0 + simple_model_plot_g1 + plot_layout(ncol=2) + plot_annotation(tag_levels="A")

ggsave(
    filename = "../images/figure-1--simple-model-plot.png",
    plot = simple_model_plot,
    width = 8,
    height = 4
)
